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Abstract 

In this work, closure of the Boltzmann-BGK moment hierarchy is accompHshed via projection 
of the distribution function / onto a space spanned by A^-order Hermite polynomials. While 
successive order approximations retain an increasing number of leading-order moments of /, the pre- 
sented procedure produces a hierarchy of (single) A^-order partial-differential equations providing 
exact analytical description of the hydrodynamics rendered by (A^-order) lattice Boltzmann-BGK 
(LBGK) simulation. Numerical analysis is performed with LBGK models and direct simulation 
Monte Carlo (DSMC) for the case of a sinusoidal shear wave (Kolmogorov flow) in a wide range 
of Weissenberg number Wi = ruk^ (i.e. Knudsen number Kn = A/c ~ \/Wi); k is the wavenum- 
ber, r the relaxation time of the system, A ~ tCs the mean-free path, and Cg the speed of sound. 
The present results elucidate the applicability of LBGK simulation under general non-equilibrium 
conditions. 
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I. INTRODUCTION 



Kinetic representations of hydrodynamics are potentially applicable to flow regimes be- 
yond the reach of classical (near-equilibrium) fluid mechanics. Nevertheless, the derivation 
and solution of high-order hydrodynamic equations for far-from-equilibrium flows with ar- 
bitrary geometry remains an open challenge. Computational methods are a valuable alter- 
native but even with the aid of efficient algorithms the solution of Boltzmann equations is a 
formidable task. Among different kinetic approaches, the lattice Boltzmann-BGK (LBGK) 
method has been able to span from scientific research to large-scale engineering applications. 
The LBGK method has two distinctive components largely responsible for its success; dis- 
cretization of velocity space and adoption of the Bhatnagar-Gross-Krook (BGK) collision 
ansatz. Decades of work have established that LBGK models correctly represent macro- 
scopic physics at the Navier-Stokes (N-S) level of approximation. On the contrary, it is not 
widely accepted in the fiuid mechanics community that high-order LBGK models provide 
hydrodynamic descriptions beyond the N-S equations. Efforts in establishing LBGK as a 
legitimate model for far-from-equilibrium flows must address two key points; the effect of 
velocity discretization errors and the validity limits of the BGK ansatz. 

The rigorous formulation of the LBGK method by Shan et al. (2006) places LBGK in the 
group of Galerkin procedures for the Boltzmann-BGK equation (BE-BGK) governing the 
evolution of the single-particle distribution /. In A^-order LBGK procedures the approximate 
solution is sought within a function space spanned by Hermite polynomials of order< A^. 
In this work, within the framework of Hermite-space approximation / G H^, we present a 
technique to systematically derive closed moment equations in the form of (single) A^-order 
partial-differential equations (PDEs). At each order of approximation, an increasing number 
of moments of / are preserved and, thus, the derived hierarchy of equations tends to the 
exact BE-BGK hydrodynamics as A^ — oo. To assess the derived hydrodynamic relations 
we perform numerical analysis with A^-order LBGK models PQ [2] and DSMC [3j for the 
case of Kolmogorov flow in a wide range of Knudsen/Weissenberg numbers (0.01 < Wi = 
t/T < 10); this free-space problem allows to remove from analysis all issues related to solid- 
fluid interaction and choice of kinetic boundary condition (e.g. diffuse scattering, bounce- 
back). Comparison of the derived equations for / G against kinetic simulations and 
previous theoretical expressions PQ H] from exact solution of BE-BGK uncovers capabilities 
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and limitations of lattice discretization and the BGK model in general non-equilibrium 
conditions. 



II. HIGH-ORDER HYDRODYNAMICS FROM BOLTZMANN BGK 

The single-particle distribution /(x, v,t) can determine all macroscopic properties (e.g. 
thermohydro dynamic quantities) observed in configuration space. In describing the fiow of 
simple fiuids we employ the velocity moments 

M(")(x,t) = y /(x,v,t)v"rfv. (1) 

The n-order moment {M*^") = ^ii i2 «„' ^fc = is a symmetric tensor of rank n and 

D is the velocity-space dimension. In similar fashion, hydrodynamic moments at local 
thermodynamic equilibrium are Mig^ = J f^'^v'^dv. The low-order moments {n < 2) relate 
to conserved quantities; namely mass, momentum, and energy: 

M(0) = Mg)= p; (2) 
M(i) = M«= pu; (3) 
trace(M(2)) = trace(Mg))= piu"^ + DO); (4) 

here we define 6 = ksT /m while T is the temperature, ks the Boltzmann constant, and m 
the molecular mass. We assume that the evolution of /(x, v, t) is governed by the BE-BGK 


df f _ f eg 

(5) 

where r is the so-called single relaxation time and the local equilibrium distribution f^'^ is 
given by 

r (^r-^^V^ 

(6) 



r''(x,v,t) = -^exp 



20 

An evolution equation for the n-order moment ([!]) can be readily obtained via moment 
integration over the BE-BGK Q: 

(^l + r^^M(") = M(:;)-rV-M("+i); n = 0, oo. (7) 

The obtained moment equation ([T]) is clearly not closed as it involves the higher-order 
moment M*^"+-'^^. 
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A. High-order hydro dynamic equations 



Leaving temporarily aside the problem of closing Eq. ([T]) let us observe that the evolution 
of M*^"^ is actually determined by all higher-order moments {M'^'^^; k > n}. From Eq. ^ 
we find that the first time derivative of M^*^^ is equal to the divergence of M'^""'"^^ i.e. the 
flux of moments one-order above. In the same way, the dynamics of M*^"+^'' is determined 
by M("+2) and so on. Climbing up the infinite moment hierarchy, one can express the 
evolution of M*^") in terms of arbitrary high-order moments {M*^""'"'^''; k > 1} after suitable 
combination of the moment equations. Multiply Eq. ([t]) by (l + t^): 

(^1 + ' M(") = + r^^ [mW - rV • M("+i)] , (8) 
and take divergence of the moment equation for the following (?2+l)-order: 

(^l + V ■ M("+i) = V ■ [M(:;+i) - rV ■ M("+2)] . (9) 

By using Eq. (joj) one can eliminate the term (l + r^) V ■ M*-""^^^ in Eq. ^ to obtain 

(l + r|)'MW = (1 + r|) - rV • M(«+i) 

+ r^V- V-M("+2). (10) 

The resulting expression, involving the evolution equations for M*^") and M'^"+^\ takes the 



form of a second-order PDE. The same procedure that lead to Eq. (10) can be applied in 
order to eliminate M^'^^^^ and iteratively performed an arbitrary number of times as the 
following higher-order moments consequently appear. After (A^ — 1) iterations we arrive to 
the general expression 



N~l , ^ . 7V-(fc+l) 

^(-^V-)Ml + r| Mir^) 

+ (-rV-)'^M("+^). (11) 
Notice here that the term (V-)^M*^""'"^) represents a tensor of rank n. The time evolution 



of the thermohydrodynamic variables corresponding to M^") is now given by Eq. (11) in 



the form of a A-order PDE. A single A-order equation of this kind implicitly involves 



the evolution of N velocity moments, i.e. those of order n to n + N — 1. Equilibrium 
moments readily computed from /^'^ (|6| are explicit function of mass, momentum, and 



energy; in solving Eq. (11) one still faces the problem of evaluating the non-equilibrium 
moment M^'^^^-' and its A^-order space derivatives. As elaborated in the next section, a 



possible way to close Eq. (11) is to express the non-equilibrium distribution / in terms of 
its leading-order moments {M*^'^-'; k < n + N} by means of finite Hermite series. 
Unidirectional shear flows. For the sake of analytical simplicity, we focus on the case of 
unidirectional shear flow u = ui with spatial gradients V = Vj = dyj and within nearly 
isothermal regime (M = u/\/6 ^1). Note that the studied unidirectional flow is exactly 
incompressible, hereinafter we adopt p = 1. The fundamental hydrodynamic variables thus 
are 

p(x,t) = 1, (12) 
u(x,t) = u{y,t)i, (13) 
e(x,t) = e + 0{M^y, (14) 

while the components of the n-order moment M*^"^ are 

^i'!l,...,*„(x,t) = j fv,,Vi^...ViJ\ =< Vi^Vi^...v,„ > . (15) 

For the studied flow the underlying distribution function must not vary along the x- and 
2;-axes {dx = dz = 0) while < Vy >=< >= 0, it follows that only the moment components 



< VxVy > {k = 0, oo) exhibit spatial variation. The iV-order equation (11) for the fluid 
velocity u{y, t) then reduces to 



r— 1 + T— u = 

dt\ dtj 

N-l ^ ^ ^ (Af-l-fc) 

,fc 



k = l ^ ^ 



< VxVy >eq 



+ (-rV)^< vxv^ >, (16) 
after recalling conservation of momentum u =< Vx >=< >eq- Hereafter, we refer to 



each A^-order PDE defined by Eq. (16) as the A^-order hydrodynamic description of the 
flow. More explicitly, Eq. (16) defines the following approximations for the studied flow: 
first-order (A^ = 1) 

_ = _V < VxVy >, (17) 



second-order (A^ = 2) 



d\du _ 
1 + r- - = -V 



+ rV' < v^v'y >, (18) 



third-order (A^ = 3) 



d\'^ du ( d 



< v^vl >, (19) 



and fourth-order (A^ = 4) 



du ( ^ d 



+ r^V^ < > . (20) 

The resulting expressions are not closed uniquely due to the presence of high-order terms 
(— rV)^< VrcVy >. If high-order terms are dominant |(tV)^| > KrV)^^"*^!, precise knowl- 
edge of the distribution / is required for accurate calculation of high-order (non-equilibrium) 



moments in Eqs. (17)-(20). On the other hand, flow regimes where |(tV)^| < KtV)^"""^! 
will permit certain approximations of / in terms of its N leading-order moments to produce 
accurate equations in closed form. 



III. HERMITE EXPANSION OF THE BOLTZMANN DISTRIBUTION 

As originally proposed by Grad (1949), the single-particle distribution can be expressed 
in terms of hydrodynamic moments via Hermite series expansion 

/(x,v,t) = r^(v)^-C(")(x,t) : H(")(v) (21) 

n=0 

with being the Gaussian weight (i.e. Maxwellian distribution for p = 1): 

^"W-(2?^«='K-^)- '''' 
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The Hermite polynomials in velocity are defined by the Rodrigues' formula: 

H(")(v) = (-l)"0teMV"e-^, (23) 
while the Hermite coefficients are 

C(")(x,t) = I /(x,v,t)HW(v)rfv. (24) 



Both H(") and C^") are n-rank symmetric tensors; the product C*^"^ : H^") in Eq. n2lh and 
hereafter represents full contraction. Each component of H^"^ (v) is an ra-degree polynomial 
in velocity v, the first four Hermite polynomials in particular are 

i/(°)(v) = 1, (25) 
H^M = ^v^, (26) 



02 

e 



H\-\^) = \{v,v,-e5,,), (27) 



and 



Hiiki^) = TT^i^^i^fc - divi6jk + Vjdik + VkSij)]. (28) 



ijk 

Hermite polynomials satisfy the orthogonality condition 



< 



jj(n) y jMjjM _ (V m ^ n) (29) 

and, hence, span the Hilbert space of square- integrable functions 5'j(v) with inner product 
< gi^Qj >= J f^Qi gjdv. Another fundamental advantage of employing the Hermite poly- 
nomial basis is that the n-order Hermite coefficient is a linear combination of the leading 
n-order moments of /. For example, 

C(o) = M(o) = p, (30) 

^5C(i) = M(i) = pu, (31) 
ec^'^) = m(2) - pel. (32) 

In similar fashion, the equilibrium distribution can be expressed as the Hermite expansion 
of the Maxwell-Boltzmann distribution (|6]) : 

oo ^ 

/^•^(x, v,t) = r^(v) ^ -CW(x,t) : HW(v). (33) 



n=0 



The Hermite coefficients Ci^^ can be readily computed using Eq. (|6|) for /^"^ in Eq. (|24[) 



A. Closure of hydrodynamic equations via Hermite expansions 



Successive order approximations can be obtained by truncating the infinite Hermite series 



(21) at increasing orders, the iV-order approximation 

^ 1 

r(x, v,t) = r^(v) J2 ^C(")(x,t) : H(")(v) 
^ — ^ n! 



n=0 



(34) 



expresses the distribution function in terms of its leading A^-order moments. The approx- 
imation / = G is tantamount to projecting the distribution function onto a finite 
Hilbert space spanned by the orthonormal basis of Hermite polynomials of order < A^. 



Due to orthogonality of the Hermite basis (29), a finite expansion (34) and the infinite series 



representation of / (21) give the same leading moments 



/v"rfv = / /^v"rfv; n<N. 



(35) 



While low order moments are preserved the higher-order moments {n > N) can be approx- 
imately expressed in terms of low-order moments. In order to close the A^-order hydrody- 



namic equations (17)-(20) we employ 



(36) 



Hence, within the framework of projection onto H^, the closed-form approximations below 
are obtained for unidirectional shear flow [see appendix |A] for detailed derivation]; / G H^: 

a \ an 



/ G e^: 



/ G H^: 



d 



d"^ \ du 



dt^ dt 



d 

1 + 3r— ] r^V^M, 
at 



(38) 



d 

1 + 3r— + 3r 



d 



dt 



+ r 



\ du 



dt^ dty dt 



(39) 



As evidenced by Eqs. (30)-(32) for {C^"^; n < 2}, second- or higher-order expansions 
(A^ > 2) are required to satisfy conservation of mass, momentum, and energy. 
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IV. iV-ORDER LATTICE BOLTZMANN BGK METHOD 



The rigorous formulation of so-called A^-order lattice Boltzmann models introduced by 
Shan et al. (2006) is based on the projection of the continuum distribution function onto 
so that /j(x, t) = /^(x, Vj,t) at a finite discrete-velocity set {v^; z = 1,Q}. Since the finite 
set of distributions = 1,<5} is expressed by A^-order Hermite series, Gauss-Hermite 

(G-H) quadrature with algebraic degree of precision d > 2N allows for exact integration 
of the leading A^-order velocity moments. Once velocity abscissae Vj and weights Wi are 
determined by a proper G-H quadrature formulae [21 E] one has 

M(")(x,t) = j /(x,v,t)v"rfv 
Q 

= J]w;,/,(x,t)v7; n = 0,N. (40) 



i=l 



Note that all Hermite coefficients (24) in the expansion of / (34) are then exactly integrated 
as well. At the same time, high-order G-H formulae determine velocity sets {vj;i = 1,Q} 
that fulfill high-order moment isotropy required for hydrodynamic representation beyond N- 
S [HI E] ■ A collateral conclusion of the Hermite expansion formulation is that the employed 
number Q of lattice velocities (i.e. quadrature points) sets an upper limit on the attainable 
order of hydrodynamic description. 

The Lattice Boltzmann-BGK Equation. The Hermite expansion formulation [2] places 
LBGK in the category of Galerkin methods, within this theoretical framework the evolution 
equations 

|^ + v,-v/^, = -^^^^ (^ = i,g) (41) 

for /j(x, t) can be systematically derived via approximation in velocity function space H^. 



The equilibrium distribution f^'^ G in Eq. (41) takes the form 

N 



/r(x,t) = /^(v.) lc£)(x,t)HW(v,). (42) 



n=0 



A. The LBGK Algorithm 



Conventional LBGK algorithms for solving Eq. (41) use an operator splitting tech- 



nique and, thus, advance in two steps: advection /f(x, t) = /j(x — VjAt,t) and colhsion 
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fii^jt + At) = (x, t) — [/f (x, t) — /f''] At/r. These steps do not constitute a standard 
Galerkin procedure, where one would directly compute the evolution of the Hermite coeffi- 
cients. As a consequence, conventional LBGK algorithms exhibit an undesired dependence 
on the flow field alignment with the underlying lattice [H HO]. This numerical anisotropy 
becomes noticeable at finite Knudsen or Weissenberg numbers where non-equilibrium ef- 
fects are important. For non-equilibrium systems will lie outside but the problem 
is effectively solved using a so-called regularization procedure [10], i.e. by re-projecting the 
non-equilibrium component /"^^ = /f — f^'^ onto H^; 

T = r (v.) ^C„e(")(x,t)H(")(vO (43) 

n=0 

where 

Q 

Cj-\^,t) = 5^^,/f(x,t)H(")(v,). (44) 



The re-projected non-equilibrium component (43) can be reintroduced at the collision step: 

/.(x + v„ t + At) = r + (i - v) 

Provided that Hermite expansions for f- (42) and /j (43) are truncated at the same 
A^th-order, the re-projection step keeps fi within as it must be the case for standard 
Galerkin procedures. The re-projection of onto is indispensable to ensure that the 
leading A^-order moments of / are exactly integrated via G-H quadrature so that simulated 
dynamics becomes independent of lattice-flow alignment. 



V. NON-NEWTONIAN KOLMOGOROV FLOW 



The decay of a sinusoidal shear wave in free space, also known as Kolmogorov flow, 
is a useful benchmark to assess derived hydrodynamic descriptions and kinetic methods 
employed in this work. In order to characterize the flow at arbitrary non-equilibrium con- 
ditions we employ the Weissenberg number Wi = t/T = rvk"^ where v = t9 is the kine- 
matic viscosity and T = vk'^ determines a characteristic decay time. Assuming a mean-free 
path A = rv^, the employed Weissenberg number directly converts to a Knudsen number 
Kn = Xk = y/Wi. In order to remain within laminar and nearly isothermal regimes the flow 
Mach number is kept small M = Uo/\/9 < 0.1; thus Re = U^/vk = Mj \/Wi < 1 is always 
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below the stability limit Re < \/2. Kinetic initial conditions are given by a distribution 
f{y, V, 0) = f'^'^ip, u{y, 0),6), i.e. local equilibrium. For this arbitrary choice of initialization 
the collision term in the kinetic equation vanishes and the simulated dynamics is collision- 
less at t = 0. As a consequence, initial conditions at hydrodynamic level are given by the 
free-molecular flow solution [T]: 



- Uo sm{ky) — exp 



2+21 



9k't 
2 



; n > 0. (46) 



We remark that after the choice of initialization at local equilibrium the microscopic dynam- 
ics remains practically coUisionless for a finite time t < r, therefore, (viscous) Newtonian 
behavior or purely exponential decay can only be observed after time intervals of the order 
of the relaxation time. The analytical description of the flow at arbitrary Wi is given by 
solution of the hydrodynamic approximations, i.e. Eqs. (37)-(39), derived in Sec. Ill via 
Hermite-space approximation / G H^. For a periodic wave, the solution to each iV-order 
hydrodynamic equation is expressed by: 

N 

u{y, t) = Y, CnIm{e'^^e-""(*+^")}. (47) 

n=l 

Each mode in the solution is determined by the complex frequencies Un(Wi) = Re{oJn} + 
ilm{un} {n = 1,A^), these values are the roots of the dispersion relation (i.e. a A^-order 
polynomial) that corresponds to the A^-order hydrodynamic approximation. The constants 
Cn and (pn in the particular solution can be determined by imposing A^ initial conditions given 



by Eq. (46) and symmetry constraints. While (positive) real roots produce exponentially 
decaying modes, each pair of complex conjugate roots describes two identical waves (i.e. 
same amplitude C and phase 0) which combine into a single standing wave that decays in 
time. 



A. Numerical simulation 



The decay of a velocity wave u{y,0) = Uq sin ky of wavenumber k = 2T[/ly is simulated 
with two different kinetic methods: the direct simulation Monte Carlo (DSMC) algorithm 



described in |i3j and the LBGK scheme described in Sec. IV A In the analysis of DSMC 
results, given that r is not a simulation parameter for this method, we use Wi ~ \vk^ jcg 
(i.e r ^ A/cs); the speed of sound c^, mean- free path A and viscosity v are determined from 
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the relations for a hard-sphere gas. For DSMC simulation we set M = 0.1 and employ a 
rather large number of particles {Np = 30000), ensembles {Nf. = 2000), and collision cells 
along ly {Nc = 500). To further reduce the statistical noise in DSMC results we perform 
spatial averaging 'w(t)/[^jj^^|^^] = / u{y,t)/u{y,0)dy over the wavelength segments ly/8- 
ly3/8 and ly5/8-ly7 /8, these quantities are presented in Figjlj For LBGK simulation we set 
M = 0.01 while the computational domain has x ly = 10 x 2500 nodes; in all cases the 
spatial resolution is conservatively larger than that determined by grid convergence tests. 
For the present results we employ the D2Q37 model (two-dimensional lattice with 37 states) 
corresponding to a G-H quadrature rule with algebraic degree of precision (i = 9 [7], i.e. 
permitting the exact integration of fourth-order moments. Different A^-order truncations of 
the Hermite expansions are implemented on the D2Q37 lattice; we refer to these schemes 
as D2Q37-H2 (iV = 2), D2Q37-H3 (A^ = 3), and D2Q37-H4 (A^ = 4). As in previous 
studies with regularized LBGK algorithms [H [10], the present results are independent of 
the flow-lattice alignment. In Fig. [T] we present the velocity field at Wi = 0.1,0.5,1,10 



given by DSMC and LBGK simulation, as well as analytical solution (47) of Eqs. (37)-(39). 
As expected, since Hermite-space approximations / G underpin the A^-order LBGK 
method, the flow simulated by LBGK models is exactly described by analytical solution to 



Eqs. (37)-(39) at arbitrary Wi. The DSMC method, which does not resort to discretization 
of velocity space nor the BGK collision ansatz, is in good agreement with LBGK and the 
/ G approximations in the parameter range < Wi < 1. 



B. Long-time decay and hydrodynamic modes 

The long-time dynamics becomes independent of the choice of initial condition for t/r = 
tuk'^/Wi 3> 1. The long-time solution of the flow is determined by the decay frequency 
u{Wi) with the smallest real part. In Newtonian regime {Wi=0), N-S solution yields a single 
hydrodynamic mode u = lm{Uoexp{iky — ut)} describing purely exponential decay with 
u = l/uk'^. Hermite-space approximations / G {N = 2,3,4) predict a long-time decay 
cu{Wi) [see Fig. |2] determined from the set of roots {cu^; = 1, A^} of dispersion relations 



corresponding to Eqs. (37)-(39). An alternative approach to Hermite-space approximations 



12 



11 

0.9 

0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
0, 



---feH'^ ' 
—feH^ - 

A D2Q37-H2 
□ D2Q37-H3 o.e 
o D2Q37-H4 
+ DSMC 




(a) 



(b) 










A"* A 

V + + ^' + + / 




+ + 


\ \ A ^ 




A V .'' ^ 


''-*.A- 




\ \i / 













(c) 



(d) 



FIG. 1: 



C/o sm(fcj/) 



VS. tiyk^; (a) VFi = 0.1 (b) Wi = 0.5 (c) Wi = 1 (d) VFi = 10. Dotted line 



(/ G H ): analytical solution of Eq. (37). Dashed line (/ G EI ): analytical solution of Eq. (38). 



Solid line (/ G M^): Analytical solution of Eq. Markers: (A) D2Q37-H2, (□) D2Q37-H3, 

(O) D2Q37-H4, (+) DSMC. 

is provided by formal solution of BE-BGK with the method of characteristics [H H]: 

/(x,v,t) = /o(x - vt, v)e"- 

+ e-7"''(x-vrs,v,t-rs)cis. (48) 
Jo 

Hydrodynamic relations for arbitrary Wi can be derived by taking velocity moments of 



Eq. (|48|); in the long-time limit t r of the studied shear flows the following dispersion 

Tu = 1 — y/n z exp(z^) erfc(z) (49) 



relation is obtained |1] 
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FIG. 2: Long-time decay: (a) ^ Markers: (A) / G [Eq. (37)], 



(□) / G [Eq. (38)], (Q) / G [Eq. (39)], (X): / G 11°° [numerical solution of Eq. (49)]. 



Dashed line: Wi ^ 1 approximation [Eq. (50)]. Solid line: Wi ^ 1 approximation [Eq. (51)] 



with z = (1—™) / y/2Wi. Numerical solution to Eq (49 ) is presented in Fig. |2| this dispersion 
relation has one trivial solution u = 1/t and a second root u = u{Wi) also on the positive 
real axis (Re{a;} > 0, Imjw} = 0). Based on asymptotic analysis of the exact solution of 
BE-BGK approximate explicit expressions have been proposed P: 

Vi + mi - 1 



and 



a; 



2Wi 



for Wi < 1, 



1 ± VI - AWi 
2Wi 



for Wi > 1. 



(50) 



(51) 



In Fig. [2| different Hermite-space approximations / G {N = 2, 3, 4) which exactly de- 
scribed LBGK results in Fig. [T] are now compared against numerical solution to the exact 



dispersion relation (49) and asymptotic approximations (50)-(51). All roots of the different 



dispersion relations have a positive real part indicating time decay of the flow, the non- 
Newtonian decay is always slower than the Newtonian decay Rejo;} < uk"^ for Wi > 
and becomes Rejcj} ~ 1/r for Wi > 1. At a first glance, the studied expressions provide 
comparable results in the limits Wi — and Wi — )■ oo while significant disagreement is 



observed for W ^ 1. Notice that Eq. (51) is the dispersion relation corresponding to the 
telegraph equation [i.e. Eq. (37)] derived for / G H^. 
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VI. CONCLUSIONS AND DISCUSSIONS 



Provided that BE-BGK is a valid model, moment equations derived for / G are in 
principle not constrained to near-equilibrium conditions. For unidirectional and isothermal 
shear flow, Hermite space approximations of different order {/ G H^; = 2,3,4} led to 
A^-order PDEs (37)-(39) for the evolution of fluid momentum [see appendix [A| for detailed 
derivation]. The studied Kolmogorov flow represents an initial value problem in free-space 
with kinetic initialization at local equilibrium, particular analytical solution to Eqs. ([37))- 
(39) has been compared against kinetic simulation via LBGK and DSMC [see Fig. [l|. We 
found that derived A^-order hydrodynamic equations predict exactly all hydrodynamic modes 



present in the flow simulated by A^-order LBGK models. We conclude that Eqs. (37)-(39) 
can be used to benchmark LBGK algorithms at arbitrary Wi and Kn number. High- 
order LBGK models and corresponding Hermite-space approximations (e.g. D2Q37-H4 and 
/ G H'^) are in good agreement with DSMC results in a wide region Wi ~ Kn^ < 1 extend- 
ing well beyond N-S hydrodynamics. These results indicate that in the region Wi < 1 the 
BE-BGK moment hierarchy approximates fairly well the low-order moments of the Boltz- 
mann equation with binary collision integral. A signiflcant disagreement exists between 
LBGK and DSMC solutions in the region Wi > 1 as seen in Fig. [l](d). 
Hereafter, we put aside a discussion on the validity of the BGK ansatz for far-from- 
equilibrium flows (e.g. Wi > 1 or Kn > 1). Instead, we proceed to study the effect of 
velocity-space discretization when solving the continuum BE-BGK over the entire parame- 



ter range < Wi < oo. The dispersion relation expressed by Eq. (49) coming from exact 
solution of BE-BGK (/ G EI°°) for t ^ r has two branches of solutions [see Fig. [2]^a-b)]. 
Meanwhile, the dispersion relation corresponding to Hermite-space approximation / G 



admit A^ roots; it follows that initial conditions may excite spurious modes in Eqs. (37)-(39). 
In order to remove initialization from analysis we examine the long-time behavior t ^ r 
characterized by the fundamental frequency u{Wi). While Re{a;} > determines the flow 
decay rate or momentum dissipation, an imaginary component lm{uj} 7^ is responsible 
for time oscillations or momentum wave propagation as observed in Fig. [TJ^c-d). We have 
compared in Fig. [2] the long-time frequency u{Wi) determined from Eqs. (j37|)-([39|) against 



u{Wi) according to Eq. (49). After truncation of the Hermite series, or corresponding ve- 
locity space discretization, dissipative properties of the flow can still be well represented for 
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Wi ^ 1, where Re{aj}/z/A;^ ~ 1, and Wi S> 1, where Re{a;}/z/A;^ ~ 1/Wi. The imagi- 
nary parts also approximate the exact BE~BGK prediction lm{u}/h'k'^ = in both hmits 
Wi — )■ and Wi — )■ oo as seen in Fig. |2](b). Notice that odd-order approximations (e.g 
/ G H^) yield a real- valued frequency u for all Wi while even-order approximations admit 
a long-time frequency with non-zero imaginary part at sufficiently high values of Wi; i.e. 
Wi > 0.25 for / G and Wi > 0.388 for / G H^. In the case of Hermite-space approxi- 
mations of even order when Wi ^ 1, time oscillations may persist in the long-time solution 
as the oscillation period becomes smaller than the decay time; e.g. Re{ci;}/Im{u;} = VWi 
for / G H^. As observed in previous work PQ |TT], a second-order approximation / G can 
be employed to model a viscoelastic response in high-frequency oscillatory flows similar to 



that observed for a Maxwell fluid and governed by the telegraph equation (37). 
LBGK methods and extensions: The LBGK method has been extensively employed for 
macroscopic description of various physical phenomena (e.g. microfiuidics, turbulence, 
reaction-diffusion, phase transition), albeit the exact (high-order) moment dynamics that 
different LBGK algorithms produce has not been fully elucidated. This inconvenience is 
partly because Chapman-Enskog (C-E) expansions, which have emerged as the preferred 
closure procedure, become increasingly difficult when carried to high-orders. The approach 
presented in this work allows to close the LBGK moment hierarchy circumventing C-E 
techniques. At the same time, it is straightforward to determine the C-E expansion order 
that correspond to a particular Hermite-space approximation [see Shan et al. (2006)]. The 



moment-equation hierarchy presented by Eq. (11) when combined with different Hermite- 
space approximations can be applied for a priori design of LBGK schemes that solve high- 
order and non-linear PDEs governing numerous complex physical systems beyond fluid me- 
chanics. It is also worth to remark that a relatively simple algorithm, based on fully-implicit 
and low-order finite-difference schemes, offering significant computational advantages can 
be effectively employed for the numerical solution of PDEs involving high-order derivatives 



in time and space, e.g. see Eq. (39) with hyperviscosity. 

The validity limits of BE-BGK: The main scope of this work is not to establish the va- 
lidity of BE-BGK in far-from-equilibrium conditions; efforts in that area could compare 
the presented analytical expressions against experimental data or more extensive numerical 
analysis via alternative methods. From results in this work it is clear that DSMC, which 
emulates the Boltzmann equation with a binary collision integral, and BE-BGK produce 
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similar solutions for the studied shear flow in the region Wi = ruk"^ < 1. Nevertheless, 
the upper applicability limit of BE-BGK for describing macroscopic physics remains to be 
established when the system dramatically departs from equilibrium conditions. 
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Appendix A: A^-Order hydrodynamic equations 

Owing to geometrical simplicity, the studied shear flow is incompressible p = 1 and 
u = u{y,t)i. 

1. Hydrodynamic Approximation in 

Approximation within space requires that all distribution functions be second-order 



Hermite expansions. Eq. (34) yields 



/ = f + luv, + A^{<vl>-e) {vl - 6) 

+ i.{<vl>-e)ivl-e) 

+ P < V^Vy > V^Vy] (Al) 

and the equilibrium distribution becomes 

r = [1 + luv^ + ^.u\vl - 9)] . (A2) 



From Eq. (A2) we obtain the equilibrium moment 

< V^Vy >eg= 0, (A3) 



while Eq. (Al) gives the third-order moment 

< v^Vy >= 9u. (A4) 
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Using Eqs. (A3)-(A4) one can close the second-order hydrodynamic description given by 

d\ du 



Eq. (18): 



This equation is known as the telegraph equation. 



2. Hydrodynamic Approximation in 



(A5) 



The / G H'^ approximation leads to 



/ = r[l + luv^ + j,{<vl>-9) {vl-d) 

+ ^{<VI>- e){vl-e) + \ < V^Vy > V^Vy 

+ w{<vl> -3ue) {vl - 3vJ) 



+ 2^3 < vlVy > {vlVy-Vy9)], 



and the equilibrium distribution 



r = + \uv, + ^.u\vl -6) + ^,u\vl - 3v, 



From Eq. (A7) one gets equilibrium moments 



< Va:Vy >eq= 0, < V^V^ >eq= Ou, 



while Eq. ( A6 ) yields the fourth-order moment 



< Vr^Vy >= 36 < VxVy > 



(A6) 



(A7) 



(A8) 



(A9) 



Recalling Eq. (17) we have V < v^Vy >= — m, and thus we can close Eq. (19): 

d"^ \ du 



d 



df^ dt 



d 

1 + 3r— ] r^V^M. 

at 



(AlO) 
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3. Hydrodynamic Approximation in 



Carrying the Hermite expansion to the fourth-order gives 



[Vr.Vy - Vxt 



f = f + luv, + ^{<vl>-e){vl-d) + ^{<vl>- e) {vl-e) + l< v^vy > v^vy 

+ (< > -'iuO) {vl - ?>v^e) + <vl> {yl- ?>Vye) + {< v^y > -u9) (? 
+ W< > i^^y - ^y(^) + 2^ (< > "6 < > +3^^') {v^ - Qvld + W^) 
+ 2^ (< < > -6 < ^1 > +3^') « - Qvp + 3^') 

+ >-<vl>9-<vl>0 + 0') [vlv'y - vie - vie + e^) 

+ el^ (< v^vl > -3 < v^Vy > e) {v^v^ - v^Vye) + {< v^Vy > -3 < v.^Vy > e) {vlvy - vi;Aj0^] 



and 



r = /""[I + luv, + i,n\vl -e) + ^u^vl - 3vJ)] + ^u\vt - 6vJ + 3e% (A12) 



Thus, Eq .(A12) yields the following equilibrium moments: 



< V^Vy >eq= 0, < V^vl > eq= 6'm, < V^V^ > eq= 0. 



(A13) 



From Eq. (All) the / G EI approximation to the fifth-order moment is 



< VxVy >= Qe < V^Vy > -35 M. 



(A14) 



Invoking Eq. (18) we have 



V* < v^v^ >-- 



6e 



l + T- 



d\' d 



dt dt 



and the fourth-order hydrodynamic description [Eq. (20)] in closed-form reads 



(9 (9^ 



dt 



\ du 



dt^ / dt 



(A15) 



dt ot^ 



(A16) 



[1] C. E. Colosqui, H. Chen, X. Shan, I. Staroselsky, and V. Yakhot, Phys. Fluids 21, 013105 
(2009). 

[2] X. Shan, X. F. Yuan, and H. Chen, J. Fluid Mech. 550, 413 (2006). 
[3] F. Alexander and A. Garcia, Comput. Phys. 11, 588 (1997). 



19 



[4] H. Chen, S. A. Orszag, and I. Staroselsky, J. Fluid Mech. 574, 495 (2007). 

[5] C. Cercignani, Mathematical methods in kinetic theory (Plenum Pub. Corp., 1969). 

[6] H. Grad, Commun. Pure Appl. Math. 2 (1949). 

[7] X. Shan and H. Chen, Int. J. Mod. Phys. C 18, 635 (2007). 

[8] H. Chen, I. Goldhirsch, and S. A. Orszag, Journal of Scientific Computing 34, 87 (2008). 

[9] H. Chen and X. Shan, Phys. D 237, 2003 (2008). 
[10] R. Zhang, X. Shan, and H. Chen, Phys. Rev. E 74 (2006). 
[11] V. Yakhot and C. E. Colosqui, J. Fluid Mech. 586, 249 (2007). 



20 



